home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 32 / Mac Magazin and MacEasy Magazine CD - Issue 32.iso / Grafik & Text / OzTeX3.0 / MetaPost / Inputs / marith.mp < prev    next >
Text File  |  1996-08-24  |  4KB  |  157 lines
  1. % Macros for dealing with very large and very small numbers in `Mlog form'.
  2.  
  3. % A number x in Mlog form represents mexp(x) if x is an even multiple of
  4. % epsilon and -mexp(x) if x is an odd multiple of epsilon, where epsilon=1/65536
  5. % is the basic unit for MetaPost's fixpoint numbers.  Such numbers can represent
  6. % values large as 3.8877e+55 or as small as 1.604e-28.  (Anything less than that
  7. % is treated as zero.)
  8.  
  9. % Mlog_str <string>            convert a string like "6.02e23" to Mlog form
  10. % Mexp_str <M_numeric>         convert from Mlog form to a string like "6.02e23" 
  11. % Meform(q)                    find (x,y) such that q is the Mlog form of x*10^y
  12. % Mlog <numeric>               convert a number to Mlog form
  13. % Mexp <M_numeric>             convert from Mlog form into ordinary numeric form
  14. % Mabs <M_numeric>             absolute value
  15. % Mlog_Str <string or numeric> convert to Mlog form; unknown results in unknown
  16. % <M_numeric> Madd <M_numeric> add
  17. % <M_numeric> Msub <M_numeric> subtract
  18. % <M_numeric> Mmul <M_numeric> multiply
  19. % <M_numeric> Mdiv <M_numeric> divide
  20. % <M_numeric> Mleq <M_numeric> the boolean operator <=
  21. % Mzero                        constant that represents zero
  22. % Mten                         constant that represents 10.0
  23.  
  24. % All other externally visible names start with `M' and end with `_'.
  25.  
  26.  
  27.  
  28. warningcheck := 0;  % Need warningcheck:=0 anytime these macros are in use
  29.  
  30. if unknown Mzero:   % If we have not already seen this file
  31.  
  32. newinternal Mzero;
  33. Mzero := -16384;    % Anything at least this small is treated as zero
  34.  
  35.  
  36. input string.mp
  37.  
  38.  
  39.  
  40. % Ideal value of mlog(10) is 589.4617838 or 38630967.46/65536.  To get an even
  41. % numerator, we round to 38630968/65536 or 589.4617920.
  42. newinternal Mten;
  43. Mten := 589.46179;
  44.  
  45. % Convert x*10^y to Mlog form, assuming x is already in Mlog form.
  46. primarydef  x M_e_ y = (x+Mten*y) enddef;
  47.  
  48. def Mgobble_ text t = enddef;
  49. pair M_iv_; M_iv_=(0,4);
  50.  
  51.  
  52. % String s is a number in scientific notatation with no leading spaces.
  53. % Convert it to Mlog form.
  54. vardef Mlog_str primary s =
  55.   save k, t, e, r;
  56.   string t;
  57.   if substring(0,1) of s="-":
  58.     r=-epsilon; t=substring (1,infinity) of s;
  59.   else:
  60.     r=0; t=s;
  61.   fi
  62.   let e = Mgobble_;
  63.   if begingroup scantokens substring M_iv_ of t endgroup >= 1000:
  64.     k = cspan(t,isdigit);
  65.     t := substring M_iv_ of t  &  "."  &  substring(4,k) of t  &
  66.          substring (k if substring(k,k+1) of t=".": +1 fi, infinity) of t;
  67.     r := r + (k-4)*Mten;
  68.   fi
  69.   let e = M_e_;
  70.   r + Mlog scantokens t
  71. enddef;
  72.  
  73.  
  74. % Convert x from Mlog form into a pair whose xpart gives a mantissa and whose
  75. % ypart gives a power of ten.  The constant 1768.38985 is slightly more than
  76. % 3Mten as is required to ensure that 3Mten-epsilon<=q+e*Mten<4Mten-epsilon.
  77. % This forces the xpart of the result to be between 1000 and 10000.
  78. vardef Meform(expr q) =
  79.   if q<=Mzero: (0,0)
  80.   else:
  81.     save e; e=floor((q-1768.38985)/Mten);
  82.     (Mexp(q-e*Mten), e)
  83.   fi
  84. enddef;
  85.  
  86. % Perform the inverse of Mlog_str, converting from Mlog form to a string.
  87. vardef Mexp_str primary x =
  88.   save p;
  89.   pair p; p=Meform(x);
  90.   decimal xpart p
  91.   if ypart p<>0: & "e" & decimal ypart p fi
  92. enddef;
  93.  
  94.  
  95. vardef Mabs primary x = x*.5*2 enddef;
  96.  
  97. % Convert a number to Mlog form
  98. vardef Mlog primary x =
  99.   if x>0: Mabs mlog x
  100.   elseif x<0: epsilon + Mabs mlog(-x)
  101.   else: Mzero
  102.   fi
  103. enddef;
  104.  
  105.  
  106. % Convert a number from Mlog form
  107. vardef Mexp primary x =
  108.   if x=Mabs x: mexp x
  109.   else: -mexp x
  110.   fi
  111. enddef;
  112.  
  113.  
  114. primarydef a Mmul b =
  115.   if (a<=Mzero) or (b<=Mzero): Mzero
  116.   else: (a+b)
  117.   fi
  118. enddef;
  119.  
  120.  
  121. primarydef a Mdiv b =
  122.   if a<=Mzero: Mzero
  123.   else: (a-b)
  124.   fi
  125. enddef;
  126.  
  127.  
  128. % 256ln(1579)=123556596.0003/65536=1885.3240356445
  129. secondarydef a Madd b =
  130.   if a>=b: (Mlog(1579 + Mexp(b Mmul (1885.32404-a))) + a-1885.32404)
  131.   else:  b Madd a
  132.   fi
  133. enddef;
  134.  
  135.  
  136. secondarydef a Msub b = a Madd (b-epsilon) enddef;
  137.  
  138.  
  139. tertiarydef a Mleq b =
  140.   (if a=Mabs a:
  141.     if b=Mabs b: a<=b else: false fi
  142.    elseif b=Mabs b: true else: b<=a fi
  143.   )
  144. enddef;
  145.  
  146.  
  147.  
  148. vardef Mlog_Str primary x =
  149.   if unknown x: whatever
  150.   elseif string x: Mlog_str x
  151.   else: Mlog x
  152.   fi
  153. enddef;
  154.  
  155.  
  156. fi   % end the if..fi that encloses this file
  157.